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Current molecular dynamic simulations of biomolecules using multiple time steps to update the 
slowingly changing force are hampered by an instability occuring at time step equal to half the period 
of the fastest vibrating mode. This has became a critical barrier preventing the long time simulation 
of biomolecular dynamics. Attemps to tame this instability by altering the slowly changing force and 
efforts to damp out this instability by Langevin dynamics do not address the fundamental cause of 
this instability. In this work, we trace the instability to the non-analytic character of the underlying 
spectrum and show that a correct splitting of the Hamiltonian, which render the spectrum analytic, 
restores stability. The resulting Hamiltonian dictates that in additional to updating the momentum 
due to the slowly changing force, one must also update the position with a modified mass. Thus 
multiple-timestepping must be done dynamically. 

I. INTRODUCTION 

The evolution of any dynamical variable W(qi,pi) is given by the Poisson bracket, 

d , s r , T-^fdWdH 8WdH\ , . 

-W {quPi ) = { W,H } = £(_ - - — -). (1.1) 

For the standard Hamiltonian, 

H ^^ = T,S^ +v{qil (L2) 

i 

the Poisson evolution equation l|l.l|l can be written as an operator equation 

dt ( rrii dqi ' dpi ) ' ^ ^ 

with formal solution 

W(t) = e^ T+v ) W(0) = [e £ ( T+ ^j V(0), (1.4) 
where T and V are first order differential operators defined by 

i i 

Their exponentiations, e eT and e eV , are then displacement operators which displace e& and pi forward in time via 

Pi 

<7i — > q% + e— and pi — > p t + eF l . (1.6) 

m l 

Each factorization of e e ( T+v ) into products of e eT , e eV (and exponentials of commutators of T and V) give rises 
to a symplectic algorithm for evolving the system forward in time. This is the fundamental Lie-Poisson theory of 
symplectic integrators, which has been studied extensively in the literatur oVi 3 . First and second order factorizationa 
of the form 

e<T+ V) ^ e eT e eV (1J) 

~ e hv e cT e hv (i. 8 ) 

give rises to the well-known symplectic Euler and the velocity- Verlet algorithm. Numerous higher order symplectic 
algorithms 3,4 - 5 - 6 ' 7,8 are also known, but only a special class of fourth order algorithms can have strictly positive time 
steps as these lower order algorithms&iS. 
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In many cases, the Hamiltonian of interest is of the form, 

V 2 

H(p,q) = j-+v 1 (q)+v 2 (q), (1.9) 

where there is a "fast" force component F\ = —d q v\ and a "slow" force component F2 — —d q V2- For example, in 
biomolecular dynamics, F\ can be the rapidly vibrating force of H-0 bonds and F2, the sum of non-bond forces. Since 
it seems reasonable to sample the slowly changing force F2 less frequently, one can factorize this Hamiltonian to first 
or second order, 

e At(T+V 1 + V 2 ) =e At(T+V 1 ) e AtV 2 ^ (1.10) 

= e iAty 2e A t (T + y l)e iA t y 2; (in) 

and solve for the fast force accurately using a smaller time step At — At/k, 

k 



e At(T+Vi) 



e ATT e ArVi 



glArVigArCT+VOglArV! 



k 



(1.12) 
(1.13) 



Thus the slow force is sampled at a multiple time steps of the fast force, At = kAr. In the context of biomolecular 
simulation, this form of the multiple-time step (MTS) symplectic algorithm was introduced by Grubmuller et aL— , 
and independently by Tuckerman et al.—. If a large time step At can be used in MTS algorithms, one can hope to 
simulate the motion of marcomolecules through some biologically significant time intervals. 

In the subsequent work of Zhou and Berne^ and Watanabe and Karplus^, this hope was dashed by the discovery 
of an intransigent instability. No matter how accurately one has solved the fast force, the MTS algorithm is unstable 
at At — tt/oji, where u>\ is the fast force's vibrational angular frequence. This has been described as a "resonance" 
instabilityASiiSiAl. However, the later numerical work of Barth and Schlick 18 clearly demonstrates that this instability 
exists at every mid-period as well, i.e., at At — wk/lhi — (n/2)T l7 where T\ is the period of the fast force, at 
n = 1, 2, 3..., and not just at n = 2, 4, 6,... Thus the notion of resonance is not a complete nor accurate description of 
this instability. In this work, we will show that this instability is fundamentally related to the non-analytic character 
of the harmonic spectrum and cannot be tamed by just multiple-timestepping the slow force. Stability can only be 
restored by a different splitting of the Hamiltonian requiring the slow force to be updated dynamically with a modified 
mass. 

In the next section, we analyze Barth and Schlick's model of MTS instability^ and show that static multiple- 
timestepping of the slow force destablizes the marginally-stable points of the fast force. In Section III, we show 
that an alternative splitting of the Hamiltonian, that of dynamic multiple-timestepping of the slow force, restores 
stability. In Section IV, we explain why the particular splitting worked in terms of the analytic character of the 
resulting spectrum. Section V generalizes MTS to the case of multiple forces. Section VI summarizes our findings 
and suggestions for large scale biomolecular simulations. 



II. STABILITY ANALYSIS OF MTS ALGORITHMS 

Barth and Schlicki^ have proposed the simplest and clearest model for understanding the MTS instability. This is 
a harmonic oscillator with two spring constants, 

Vx{q) = ^X iq 2 , v 2 (q) = l -\2q 2 - (2.1) 

Their numerical work unambiguously demonstrated the existence of MTS instability, but they did not carry their 
analysis far enough to pinpoint its origin. We will first complete their analysis of the symplectic Euler MTS algorithm. 

Each operator e ArT , e ArVl , when acting on the canonical doublet (p,q), produces a symplectic transformation, or 
map, 

p n+l 
r.n+1 



= e^ * =V(Ai,Ar) ' , (2.2) 



C-l )=e ArT ( P ' )=T(m,Ar)( P ' n ), (2.3) 



q T I \ q I v \ 1 
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where T and V are matrices given by 

T(m, At) 
V(A,At) 

The Jacobian of the transformation defined by 

M = 



1 

Ar/m 1 

1 -AtA 
1 



d(p n+ \q n+1 ) 
d(p n ,q n ) 



satisfies the defining symplectic condition 

M T JM = J, with J = 



-1 

1 



ensuring that detM T detAi^l. For a sequence of symplectic maps, by the chain-rule, the Jacobian multiplies 

d{p n ,q n ) _ d(p ni g n ) d{p 2 ,g2) d(pi,<?i) 
d{pa,qo) d{p n -i,q n -i)'" d{pi,qi) d{p ,q )' 



(2.4) 
(2.5) 

(2.6) 

(2.7) 

(2.8) 



Regarding 12.212.3(1 as numerical algorithms, the Jacobian matrix is just the error amplification matrix. However, 
only in the present case of linear maps I|2.2I2.3[) is the Jacobian the same as the transformation matrix itself. 
The error amplification matrix corresponding to the symplectic Euler MTS algorithm 



^At(T+Vi+V 2 ) 



e ArT e ArVi 



0(At 2 ) 



is therefore (corresponding to Barth and Schlick's Ai), 



At At n k 

T(m, T )V(A!,— )] V(A 2 ,At). 



(2.9) 



(2.10) 



The symplectic matrices T and V as defined by l|2.4|l and 1)2. 5fl . can also be expressed as exponentials of traceless 
matrices: 



T(m, At) = exp 
V(A, At) = cxp 



At 
At 




1/m 

-A 




(2.11) 
(2.12) 



For large multiple k, the fast force term in (|2.1U|) can be evaluated analytically. Using the exponential forms for T 
and V, and invoking Trotter's theorem, 



lim exp 



At 
k 



exp 




1/m 

cos(wi At) 



AtfO -Ai 




= exp 



At 



-Ai 
1/m 



mui sin(uiAt) \ _ 



(muoi) 1 sin(w 1 At) cos(oj 1 At) J 



= H(m,o;i,At), 



(2.13) 



where we have defined the fast force angular frequence uj\ = \J Ai/m. Note that one starts with Ai and m, but the 
dynamics of the system is governed by the square root oj\. Since oj\ is a non-analytic function of Ai and m, it can 
only be extracted in the limit of k — *■ oo. 

The eigenvalues of the fast force error matrix l|2.13|l is given by 



ei,2 



C ± Vc 2 - 1 



(2.14) 



with C = cos(a;) and x — uiiAt. The algorithm is marginally stable at all time step At with (ei^l = 1, but closest to 
being unstable at x = n7r, where the two eigenvalues are degenerate, purely real, and equal to ±1. 
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The error matrix corresponding to Euler MTS algorithm (|2.10|l is therefore 

e J s = H(m,wi,Ai)V(A a ,At). (2.15) 
The eigenvalues are still given by H2.14[l . but now with C altered to 

C = cos(a;) — —ax sin(x), (2-16) 

= A(x)cos{x + 5(x)), (2.17) 

with a = A2/A1, amplitude 

A{x) = y/1 + (ax/2) 2 , (2.18) 

and phase shift S(x) = tan _1 (aa;/2). The two C-functions, together with the amplitude functions ±A(x), are plotted 
in Fig.l. The Euler MTS algorithm is unstable whenever |C(x)| > 1. As shown in Fig.l, the effect of A2, no matter 
how small, is to destablize marginally stable points of the fast force Ai at x = mr into a finite band. The first band 
at x = 7T, is the half period barrier. The bands are very narrow if A2 << Ai. Within these instability bands, the 
extremes of the eigenvalues at x + S(x) = mr, when C = ±A(x), are given by l|2.14ll . 

e(x) = ± (ax/2 + y/l + (ax/2) 2 ) . (2.19) 

This is the linearly growing envelope of eigenvalues observed numerically by Barth and Schlick— . Since the eigenvalue 
departs from unity linearly as a function of x, we can characterize this instability as first order in x. This is the most 
important characterization of MTS algorithms and is plotted in Fig. 2. As one can see, as long as a is not zero, 
the departure from unity will be significant at x = %, which explains the persistence of the half period barrier. We 
emphasize that e(x) only gives the correct eigenvalues at x + S(x) = mr, when C — ±A(x). For a << 1, this means 
that e(x) is only correct at x sa mr. For other values of x, e(x) is not the correct eigenvalue and the algorithm is 
actually stable. 

The error matrix for the second order Verlet-like MTS algorithm, 

ev = V(A 2 , i At) H(m, w x , At)V(A a , ^At) (2.20) 

has the same C-function (|2.16() and therefore the identical first order instability problem. This is a surprise. As we 
will see later in Section IV, increasing the order of static MTS algorithms does little to increase its stability. 



III. RESTORING STABILITY VIA DYNAMICAL MTS 



The MTS algorithm in the last section splits the Hamiltonian as 

HM={^\^)+\^\ (3.1) 

where the parenthesis describes the full dynamics of spring Ai. This leaves A2 as only a static force with no direct 
role in changing the particle's position. We shall refer to this as static multiple-timestepping. This is not an equitable 
splitting, nor the only one possible. The Hamiltonian can alternatively be splitted as 

"<"» = (2t + 5 A " 3 ) + (^ + ^ 2 )' (3 - 2 ' 

with the constraint 

-L + -L.1. (3.3) 
mi TO2 m 

Now both springs are fully dynamical and we can use the freedom in the choice of mi and m 2 to maximize stability. 
We shall refer to this as dynamic multiple-timestepping. The Euler splitting algorithm of l|3.2|) in operator form is 

e At{T 1 + V 1 +T 2 + V 2 ) ^ e At(T 1 + V 1 ) e At(T 2 + V 2 ) (3 4) 



Consider first when both are evaluated exactly as in l|2.13[l . then the error matrix is 

e DE = H(mi, fix, At)H(ma, fi 2 , At), (3.5) 

with 



ill = \ — - an d ^2 = \ — — • (3-6) 
V wi V m-2 



The corresponding C-function is 



C = cosffftx + n 2 )At) - ( TO i^i rn 2 n 2 ) 2 g j n /Q i ^\ sin (fi 2At ). (3.7) 

The destablizing sine function term can be eliminated by choosing 

77iif2i = m,2f22 - ► m\\\ = mi\i. (3-8) 
Thus stability can be fully restored in this splitting with the choice of 

— = t ~~\ and — = — — . 3.9) 

mi Ax + A2 m m 2 X\ + A2 m 

For this choice of mi and 1112, we observe that 

Ol = - Al — fi and fl 2 = . A2 , Q (3.10) 
Ai + A2 Ai + A2 

where 



n = J±±^ (3.11) 

V m 

is the exact angular frequence of the system. This means, however that 

f7 = f2i + fi 2 , (3.12) 

i.e., the choice of mi and m 2 which restores stability also linearizes the spectrum. To compare with the static case, 
we also note that 



w Ai+A 2 m VTToT 
and 

fi 2 =afix. (3.14) 

For MTS algorithms, we do not want to evaluate the second spring force exactly, but only sparingly. Thus we 
further approximate ({3.4(1 by 

e At(Ti+Vi+T 2 + V2) ^ e At(Ti+Vi) e AiT 2e AiV 2 /g ^\ 

This is the dynamical Euler MTS algorithm with error matrix 

e DE = H(mx, n x , At)T(m 2 , At)V(A 2 , At). (3.16) 

The resulting C-function is 

C = cos(a;')(l - ^(qi') 2 ) - ax' sin(x'), (3.17) 

where 

x' = QiAt = x/Vl + a and ax' = tt 2 At. (3.18) 
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This C-function is cos(f2iAi + fl 2 At) correct to second order in il 2 At. The corresponding amplitude and eigenvalue 
functions are 



A(x') = y/1 + (ax') 4 /4, 



(3.19) 



e(x') = ± (ax') 2 /2 + a/I + (ax') 4 / 4 



(3.20) 



Thus by allowing A2 to be dynamical, the same effort in force evaluation improves the instability to second order. 
This is shown in Fig. 2. However, one can do even better. By (|3.5|l . the algorithm's stability will continue to improve 
with improvements in solving A2's dynamics. With still only one slow force evaluation, one can solve A2's dynamic to 
second order with error matrix 



C-function 

amplitude 
eigenvalue 



e DE2 = H(mi, Qi, At)T(m 2l ^At)V(A 2 , Ai)T(m 2 , ^At), 



C = cos(a;')(l - \{ax'f) - ax' {I - \{ax'f) sin(x'), 
2 8 



A(x') = y/1 + (ax'/2) 6 , 



e(x') = ± (ax'/2) 3 + ^l + (ax'/2f 



(3.21) 
(3.22) 
(3.23) 
(3.24) 



and improve stability to third order! In sharp contrast to the static case, where the use of a second order algorithm 
for the slow force yielded no improvement in stablity, the improvement here is dramatic. As shown in Fig. 2, even for 
a as large as 1/20, this second order algorithm is basically stable at x — tt. 

If one is willing to evaluate the slow force more than once, further systematic improvments are possible. The second 
spring's motion can be solve to fourth order using forward symplectic algorithm 4 A?' 1 ? : 



,At(T 2 +V 2 ) 



C2 



Aty 2e iAtT 2o |A t V 2a iAtT 2o iAty 2 



'e 3 "e 2 



e e - 



0(Atf 



Here V 2 = V 2 + -joAt 2 [V 2 , [T 2} V 2 ]\. The double commutator modifies the original spring constant A 2 to 



A 2 = A 2 (l - ±-^At 2 ) = A 2 (l - ±-{ax'f) 



24 m 2 



24 



The resulting error matrix is 



e 4A = H(mi, fii, At)V(A 2 , ^Ai)T(m 2 , ^Ai)V(A 2 , ^Ai)T(m 2 , ^Ai)V(A 2 , ^At), 

6 2 3 2 6 



(3.25) 



(3.26) 



(3.27) 



with C-function 



C 



cos(x') 



ax'-^ax'^ + ^^x') 5 



1 



'\7 



10368 



(ax') 



sin(x'), 



amplitude 



and eigenvalue function, 



A(x') 



1 + p(a*72)i° - |(^'/2) 12 + ^(ax'/2) 14 , 



z{x') = ± [V^V) - 1 + A{x') 



(3.28) 



(3.29) 



(3.30) 



The instability is now pushed back to fifth order in x. Fig. 2 shows that even for a as large as 1/20, this algorithm is 
now basically stable out to x w 677. For a = 1/400, as considered by Barth and Schlick, this algorithm has e < 1.00001 
at x ~ 50tt. There is no doubt that one has overcame the half-period barrier at x = tt. 
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IV. STABILITY EXPLAINED 

The poor stability of static multiple-timestepping can be traced to the non-analytic character the spectrum. The 
system's exact angular frequence is 



Ai + A 2 



n = y ~ - = yful + u% = viVT+^, (4.1) 

with exact C-function 

C = cos(flAt). (4.2) 



In terms of x — u>\At — y Xi/mAt and a = A2/A1, this function is non-analytic in a, 



C = cos(x\/l + a). (4.3) 

When expanded in terms of a, it has the form 



C = cos(x) — — xsa\(x)a + 



-xsin(a;) — —x 2 cos(x) 
8 8 



of + ■ ■ ■ . (4.4) 



The first order term is precisely the first order result H2.16[l . If one were able to reproduce this expansion, one could 
in principle systematically restore stability. Unfortunately one cannot; when regarding A2 as static, one must expand 
in powers of V(A2, At) oc A 2 At oc ax, and can never reproduce the term oc xa 2 in l|4.4[l in any finite order. Worse, 
second and fourth order algorithms do not even reproduce the (xa) 2 term with the correct coefficient. 
By contrast, in dynamical multiple-timestepping, one has, 

n = fL 1 +Q, 2 , (4.5) 

and the spectrum is linear in O2 - The corresponding C-function 

C = cos(OiAt + n 2 At) = cos(x' + ax'), (4.6) 

as shown in the last section, can be systematically reproduced order by order in (ax 1 ). Thus dynamical multiple- 
timestepping linearizes the spectrum and can overcome the half period barrier by going to higher order. 

V. GENERALIZATION TO MANY FORCES 

For more than two forces, the generalization is easy. Again, using the harmonic oscillator as an illustration, the 
"iV-forces" case of 

2 1 N 

i=l 



can be dynamically splitted as 



with the primary constraint 



N / 2 1 \ 



N 1 1 



£-"■ (5-3) 

and the pair- wise stability conditions, i =/= j, 



mi m 

2=1 



rriiXi = rrijXj. (5-4) 
Both can be easily satisfied by the following generalization of (|3.9|) . 

V 1 (5.5) 



rrii y]j—i Aj m ^ ^ ' m 

Thus the inverse of the dynamical mass should be chosen in proportional to the strength of the force, or the square 
of its angular frequence. 
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VI. CONCLUSIONS 



In this work, we have given a detailed analysis of Barth and Schlick's model of MTS instability 1 ^. We show that 
the instability of static MTS algorithms can ultimately be traced to the non-analytic character of the underlying 
spectrum. Static MTS algorithms are simply very poor starting points for solving such a spectrum, even if one were 
to modify or average over the slow forced. By contrast, dynamic MTS algorithms linearize the spectrum, render it 
analytic, and can improve stability systematically order by order. The use of Langevin dynamics to damp out the 
instability 20 simply masks the true dynamics of the system without fundamentally solving the instability problem. 

Realistic biomolecular simulations are too complicated for a detailed stability analysis as in the harmonic oscillator 
case. Nevertheless, the harmonic oscillator captures the essence of its fast vibrating modes. Thus the insight of 
dynamic multiple-timestepping can be applied easily. The key idea is to decompose 



m mi rri2 

and update particles affected by the slow force dynamically with mass m^. In the harmonic oscillator case, 1/mi and 
l/m,2 are to be determined in proportional to the strength, or the square of the frequence, of the force. For realistic 
simulations, one can simply determine the optimal by trial-and-error subject to the constraint i|6.1[l . 
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Figures 




FIG. 1: The fundamental C-function for determing the stability of Multiple-Timestepping (MTS) algorithms. The dashed line 
gives the stable C-function for the fast force alone, C(x) = cos(a:) where x = u>iAt, and uj\ is the vibrational angular frequence 
of the fast force. The solid lines give the C-function for the static Euler MTS algorithm, Eq. 12.16ll . To make the unstable 
regions visible, a large value of a = A2/A1 = 1/10 is used, where Ai and A2 are the force constant of the fast and slow force 
respectively. The algorithm is unstable whenever |C(x)| > 1. The most unstable point in each unstable band near x ~ nn 
touches the amplitude envelope ±A(x), Eq. l2.18H . 
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FIG. 2: The magnitude of the error matrix's eigenvalue for various MTS algorithms. The dashed line is the static Euler MTS 
algorithm. The three solids lines are the three dynamic MTS algorithms described in the text. The algorithm is unstable 
whenever \e(x)\ > 1, however, in this graph, only values at x + 8(x) = nir are true eigenvalues. See text for details. A large 
value of a = 1/20 is used to make the fourth order dynamic MTS result visible. 



